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ABSTRACT 



In many brain diseases it can be qualitatively observed that spatial patterns in blood oxygenation level dependent 
(BOLD) activation maps appear more (diffusively) distributed than in healthy controls. However, measures that 
can quantitatively characterize this spatial distributiveness in individual subjects are lacking. In this study, we 
propose a number of spatial heterogeneity measures to characterize brain activation maps. The proposed 
methods focus on different aspects of heterogeneity, including the shape (compactness), complexity in the 
distribution of activated regions (fractal dimension and co-occurrence matrix), and gappiness between activated 
regions (lacunarity). To this end, functional MRI derived activation maps of a language and a motor task were ob- 
tained in language impaired children with (Rolandic) epilepsy and compared to age-matched healthy controls. 
Group analysis of the activation maps revealed no significant differences between patients and controls for 
both tasks. However, for the language task the activation maps in patients appeared more heterogeneous than 
in controls. Lacunarity was the best measure to discriminate activation patterns of patients from controls (sensi- 
tivity 74%, specificity 70%) and illustrates the increased irregularity of gaps between activated regions in patients. 
The combination of heterogeneity measures and a support vector machine approach yielded further increase in sen- 
sitivity and specificity to 78% and 80%, respectively. This illustrates that activation distributions in impaired brains 
can be complex and more heterogeneous than in normal brains and cannot be captured fully by a single quantity. 
In conclusion, heterogeneity analysis has potential to robustly characterize the increased distributiveness of brain 
activation in individual patients. 

© 2014 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license 

(http://creativecommons.Org/licenses/by-nc-nd/3.0/). 



1. Introduction 

Brain activation maps derived from functional MRI measurements 
are usually compared between diseased and control subject groups to 
draw statistical inferences on aberrant activation patterns related to 
the neurological disease conditions. Such inferences on aberrant brain 
activation patterns rely on the spatial overlap of activated regions 
among the members of a well-defined sub-population (Haxby et al., 
2001 ; Manoach et al., 2000). However, even in the normal brain, cog- 
nitive functions depend on network activity of which the activation 
patterns are inherently distributed and are likely affected by natural 
variability of the brain's functional organization, especially in more 
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complex tasks. Any abnormalities in brain activation may not only 
express as aberrant activation levels but also as modifications in the dis- 
tribution of activated regions. Inferences on abnormal brain activation are 
therefore not unambiguous as a priori we do not know whether differ- 
ences in group averaged activation results are due to differences in activa- 
tion level or reductions in overlap of activated regions. To provide further 
insight in aberrant activation maps of impaired brains, novel measures are 
required that characterize, and preferably quantify, the heterogeneity of 
activation patterns at the individual subject level. Such measures should 
ideally characterize the organization of the distribution of brain activation 
rather than the local degree of activation overlap. 

Multiple measures for the characterization of spatial heterogeneity 
have previously been explored in other fields of medical imaging 
(Bright et al., 2009; Damona et al., 2008; Huettel et al., 2004). In the 
field of tumor imaging, spatial heterogeneity of contrast enhanced 
structural scans has for a long time been considered as a marker of ma- 
lignancy (Alic et al., 2006; Alic et al., 201 1 ; Jackson et al., 2007; Mohajen 
et al., 2011; O'Connor et al., 2011; Rose et al., 2009; Tixier et al., 2011). 
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Within the domain of fMRI, heterogeneity in time-series data has 
also been explored previously. Zang et al. (2004) presented a method 
for fMRI data analysis based on regional homogeneity, in which 
Kendall's coefficient of concordance (Baumgartner et al., 1999) was 
used to compare the similarity of time-series in a voxel to the series 
from its neighbors. This method has also been applied in studies of 
Alzheimer's disease (Liu et al., 2008), neuromyelitis optica (Liang 
et al., 2011), aging (Wu et al., 2007), individual intelligence (Wanga 
et al., 201 1 ), autism spectrum disorder (Shuklaa et al., 201 0), depression 
(Yao et al., 2009), and attention-deficit hyperactivity disorder (ADHD). 
Leech and Leech (2011) quantified spatial heterogeneity in fMRI as a 
variation in the intensities of activated and neighboring voxels. These 
publications suggest that there could be increased activation map het- 
erogeneity in certain patients. All the studies have investigated hetero- 
geneity in activation patterns between subject groups, whereas this 
study aims to quantify spatial heterogeneity on the subject level. 

Previous studies have tried to quantify the existence of heterogeneity 
in general. The current study shows how to characterize heterogeneity 
more specifically into different dimensions. We have evaluated well- 
known heterogeneity measures of functional MRI derived activation pat- 
terns in patients with epilepsy and language impairment and compared 
these to those obtained in age-matched healthy controls. The measures 
are applied to activation maps pertaining to language and motor tasks. 
The primary hypothesis is that the spatial heterogeneity of brain ac- 
tivation maps is stronger for patients than for controls. To further 
strengthen this hypothesis in relation to the disease characteristics, 
we secondarily investigated whether the increase in spatial hetero- 
geneity is specific to the language-related, and not motor-related, 
activation pattern in (language-impaired) patients. The focus of 
the current study is on the application and evaluation of heteroge- 
neity measures that can be applied to activation maps of individual 
subjects with a neurological disorder. 



2. Methods 

2.1. Overview 

The selected heterogeneity measures focus on different aspects 
of spatial distributiveness, including the shape (compactness), com- 
plexity in the distribution of activated regions (fractal dimension 
and co-occurrence matrix), and gappiness between activated 
regions (lacunarity). The various heterogeneity measures are de- 
scribed below and are calculated for the activation maps of a lan- 
guage and a motor task. Furthermore, a support vector machine 
(SVM) classifier is used to analyze the discriminative power of the 
heterogeneity measures. The block diagram in Fig. 1 outlines the 
methodological steps used for the spatial heterogeneity analysis. 

To exclude the influence of other factors, image noise and head mo- 
tion analysis is applied to the activation maps. Head motion analysis is 
performed by computing displacement, entropy and smoothness of 
the three dimensional head movement time-series. The noise level is 
determined by the signal-to-noise-ratio (SNR) from the acquired 
echo-planar images. 

2.2. Heterogeneity analysis 

In this section we explain each method to determine the different 
heterogeneity measures in detail. In Appendix A the heterogeneity mea- 
sures are illustrated using a conceptual set of images. 

2.2.J. Overlap of activation 

To obtain a measure for the degree of spatial activation overlap, we 
calculated the number of overlapping (i.e. commonly) activated voxels 
in the patient and control groups separately as a function of activation 
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Fig. 1. Block diagram of the subsequent methodological steps in the heterogeneity analysis. 
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t- threshold. Note that overlap of activation represents a measure of het- 
erogeneity at the group level rather than at the individual level. 

2.2.2. Number of activated regions 

The number of activated regions is compared between patients and 
controls as a function of the r-threshold, which controls the most signif- 
icantly activated x-% (x varies from 3-10) of voxels. 

2.2.3. Fractal dimension 

A fractal dimension is a statistical index of complexity comparing how 
details in a pattern change with the spatial scale at which it is measured 
(Dougherty and Henebry, 2001). The concept of fractal dimension was 
first proposed by Mandelbrot (1983) to describe the geometry of objects 
in nature. The fractal dimension can be seen as a number that describe the 
texture (and thus heterogeneity) in the object. Objects with high fractal 
dimension are heterogeneous. There are various types of fractal dimen- 
sion analysis techniques that could be used to describe the complexity 
(fractal pattern) of objects such as the frequently used box dimension, 
also called the Renyi dimension (Mandelbrot, 1983). The box dimension 
is a computationally simple method to quantify heterogeneity in binary 
images. It is computed by imposing regular grids of a range of spatial 
scales on the object studied and counting the number of grid elements 
(boxes) n that are occupied by the object. By plotting the number of occu- 
pied boxes as a function of the reciprocal scale on a log-log axis the box 
dimension represents the slope of the regression line. 

2.2.4. Shape of activated regions 

The shape of the activated regions is compared to a ( single) perfectly 
spherical region. The measure that quantifies similarity of a region to a 
sphere is the isoperimetric quotient 0_ and is defined as 



n _ m 13 v 

in which V and S are the volume and surface area of the activated re- 
gions, respectively. The volume is the total number of all activated 
voxels and the surface area is the number of activated voxels which 
are located at the surface of the activated regions. Voxels are classified 
to be on the surface of activated regions, if at least 1 out of their 26 
neighboring voxels is not activated. The value of Q decreases when the 
region deviates from a sphere, because of the relative increase in surface 
area. Also the splitting of one large region into multiple smaller (distrib- 
uted) regions, with preservation of total activation volume, will lead to a 
lower (averaged) Q. Note that the spatial distribution of multiple (non- 
overlapping) activated regions has no effect on Q. 

2.2.5. Co-occurrence matrices 

The co-occurrence matrix is defined as the distribution of 
co-occurring intensity values (r-values) at a given spatial offset 
(Amadasun and King, 1989; Thibault et al., 2009). This matrix is widely 
used for the characterization of texture in gray-scale images (Thibault 
et al., 2009). It is computed as 

r a ii-T-T-V J 1 '^ I{P,Q^) = iandl(p + Ax.q + Ay.r + Az)=j 
<-m,a,M''J> - ^ ^ £j ^ 0, otherwise. 



Here C is the co-occurrence matrix, n, m and o represent the sizes of 
the image and /(p,q,r) is the value at voxel position (p,q,r). The co- 
occurrence matrices used in this study are computed by taking mean of 
co-occurrence matrices computed at different offsets (Ax = 1, Ay = 0, 
Az = 0), (Ax = 0, Ay = 1, Az = 0) (Ax = D, Ay = D, Az = D), 



where D is selected such that difference in entropy at D and D— 1 is less 
than 0.0001. Two texture measures, i.e. homogeneity and entropy (Filho 
and Sobreira, 2008), were derived from these co-occurrence matrices. 
Homogeneity is the measure of image uniformity, which is inversely 
related to the heterogeneity, and is defined as 

Homogeneity = T] „ < "^,!'^. 

tfl + |i-j|- 

C(y) is the value of the co-occurrence matrix at the ith row and jth 
column. Entropy is a measure of randomness in an image and increases 
with heterogeneity and is defined as 

Entropy = -EC (i,j)ln(C'(i,j)). 

ij 



Here C'(y) is the normalized value of the co-occurrence matrix at (ij) 
by dividing each element of the co-occurrence matrix by the sum of all 
matrix elements C(ij). 

2.2.6. Lacunarity 

Lacunarity is a measure of how sparse patterns fill space. It gives 
higher values where patterns have more, larger or more variable gaps 
(Dougherty and Henebry, 2001; Filho and Sobreira, 2008). Lacunarity 
was also originally developed by Mandelbrot ( 1 983 ) to describe a prop- 
erty of fractals, and has since been extended to describe real data sets 
that may or may not resemble fractal distributions. Mandelbrot proved 
that different textures might have the same fractal values, but different 
lacunarity values. Therefore lacunarity can be interpreted as a comple- 
mentary measure to the fractal dimension. Fractals depend on the geo- 
metric structure of an object whereas lacunarity is more sensitive to the 
sizes and distribution of gaps between the geometric structures. Geo- 
metric objects with low lacunarity are homogeneous because all gap 
sizes are the same, whereas high lacunarity objects are heterogeneous. 
It is important to note that objects that are homogeneous at a particular 
scale might be more heterogeneous at other scales, therefore lacunarity 
may depend on the spatial scale of measurement. 

Lacunarity can be implemented in 3D using the sliding box 
approach. According to this algorithm a box of size r x r x r slides ( over- 
lapping) over a volume image. The mass M of an image is then defined 
by the number of activated voxels in a gliding box. n(M,r) represents the 
number of gliding-boxes with size r x r x r containing M activated 
voxels. The probability distribution Q(M,r) is obtained by dividing 
n(M,r) by the total number of boxes B: 



Lacunarity at scale r is defined as the mean-square deviation (i.e. 
second moment) of the variation of mass distribution probability 0_(M, 
r) divided by its square mean: 

J2 M 2 Q_(M,r) 



£ MQ(M. r) 

M 



Conceptually, this resembles the (relative) variability of activated re- 
gion sizes normalized to the (squared) average region size at a certain 
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spatial scale r. We have used multiple cubic gliding box sizes ranging 
from 8x8x8 mm to 16 x 16 x 16 mm. 

2.2.7. Intelligent parameter selection and classification 

Heterogeneity is a complex feature that may not be captured by just 
one measure but needs multiple measures. To identify and combine the 
most relevant heterogeneity measures, a classification method was 
applied. Following Chen and Lin (2005), the F-score has been used to 
compute the importance of heterogeneity measures for distinguishing 
patients and controls. It is defined as 

F , ft (x P l-\) 2 + (x C l-\) 2 

' = ~1 P 2 1 C c 2 

p—^^2k=l( x k,i~ X i) ^k=l( x k,i~ x i) 



Here, x k , k = 1 T, is a training vector and T is total number of sam- 
ples. P and C are the total number of patients and controls, respectively. 
Xi , X? , x\ are the mean of the ith measure of the total, patient, and control 
data sets, respectively; x\ ( , is the ith measure of the fcth patient, and x c k ( , 
is the ith measure of the fcth control. The numerator indicates the dis- 
crimination between the patients and controls (interclass variation), 
and the denominator indicates the discrimination within each of 
the two sets (intraclass variation). Higher F-scores indicate better 
discrimination. 

Given a set of parameters for each subject, our objective is to classify 
each subject as a patient or control based on the values of the measures. 
To perform the classification based on the calculated heterogeneity 
measures, we have used the support vector machine (SVM) classifica- 
tion algorithm (Kumar, 2004; Zeng et al., 2007). SVM has been applied 
with a polynomial kernel (degree 2), using the selected measures 
with the F-score method. A significant advantage of the SVM algorithm 
is the search for a global and unique solution. 

2.2.8. Statistical analysis 

Differences in heterogeneity measures between patients and 
controls were assessed using a non-parametric Wilcoxon rank test 
(Conover, 1 980). A non-parametric test has been used because distribu- 
tion of the data is unknown. To determine the degree to which extent 
the different heterogeneity measures were independent of each other, 
the (Pearson) correlation matrix over all measures was calculated. 

Furthermore, the sensitivity and specificity of the different heteroge- 
neity measures were calculated to determine the potential of the differ- 
ent measures to distinguish the activation maps of patients from 
controls. This analysis was performed using 5-fold cross validation 
(Kohavi, 1995). For each fold, all the data were randomly assigned to 

five sets SI, S2 and S5 so that all sets were of equal size and each 

set contains an equal number of patients and controls. Then training 
was performed on SI, S2, S3, and S4 and testing on S5, followed by 
other combinations such that each set has been used once for testing. 
This approach has the advantage that the training and testing sets 
were both relatively large, and each data set is used for both training 
and validation on each fold. 

Receiver-operator-characteristic (ROC) analysis of all heterogeneity 
measures was used to determine the discriminative power of each 
heterogeneity measure. For each heterogeneity measure, the area- 
under-curve (AUC) was calculated by varying the thresholds. The AUC 
was also obtained using 5-fold cross validation to provide the uncertainty. 

The spatial correspondence between the activation maps of the two 
groups was characterized by the Jaccard index, which is the measure of 
similarity. This index was calculated as a ratio of the number of overlap- 
ping activated voxels to the total number of distinct voxels. Jaccard 
index will be lower for the group with less overlapping voxels, i.e. 
group with more spatial variability in activation. 



2.2.9. Motion artifacts 

To test whether patients give rise to stronger motion induced con- 
founds in functional images, for instance spurious activation regions, 
the absolute and relative head motion time-curves were evaluated in 
terms of displacement, smoothness and entropy and correlated to the 
activation heterogeneity measures. The description and results of this 
analysis are given in Appendix B. 

22.10. Noise analysis 

SNR was determined on the original T2*-weighted images for each 
subject. The region inside the brain is considered as the signal and 
outside the brain as the noise. SNR is computed by dividing mean 
value signal inside the brain by the standard deviation of the signal out- 
side the brain. 

2.3. Subjects 

Twenty-three children were recruited from the database of a 
specialized tertiary referral center for epilepsy upon having a clinical di- 
agnosis of Rolandic epilepsy. Rolandic epilepsy is the most common 
benign childhood epilepsy, has a genetic basis and is characterized 
by centrotemporal spikes on the electroencephalogram (EEC) 
(Panayiotopoulos et al., 2008). Seizures occur mostly during the 
night and involve hemifacial spasms and speech arrest. In recent 
years evidence has accumulated that Rolandic epilepsy is associated 
with co-morbidities such as language impairment, despite its mild 
seizure semiology. 

The average age at testing was 1 1.4 years (range 8-14 years) and we 
included 7 girls and 16 boys. A control group of 21 age and gender 
matched healthy children was also included (10 girls, 11 boys, mean 
age 10.3 years, range 8-14years). 

A previous study on these groups revealed a bilateral activation 
pattern with no clear (overlapping) differences in activation maps be- 
tween the children with Rolandic epilepsy and healthy peers (Besseling 
et al., 2013). 

2.3. J. Language assessment 

The Clinical Evaluation of Language Fundamentals (CELF-4) test for 
children, Dutch edition (Paslawski, 2005; Semel et al., 2010), was used 
to assess language performance (Overvliet et al., 2013). The central 
outcome measure of the CELF test is the core language score, which 
can serve as a screening measure for language impairment. The 
patients' core language score was 92 ± 18(mean ± SD), which is signif- 
icantly lower than the norm value of 100 ± 15(p = 0.047). 

2.4. MR1 

For fMRI a blood oxygenation level dependent (BOLD) sequence 
(T2*-weighted) was used at 3 Tesla field strength with the following 
settings: single-shot echo planar imaging (EPI), echo time/repetition 
time (TE/TR) 35/2000 ms, pixel size 2x2 mm 2 and 4-mm thick axial 
slices. Each scan consisted of 195 whole-cerebrum dynamic acquisitions. 
The acquisition time was 6.5 min. 

2.4.1. Language paradigm 

A standard block design was used consisting of 6 task condition 
blocks interleaved with baseline condition blocks. Each block lasted 
30 s and each paradigm started and ended with a baseline block. 

The task comprised a word generation paradigm. Subjects had to 
covertly generate as many words as possible starting with a visually 
presented letter (U-N-K-A-E-P). The paradigm consisted of 6 word 
generation blocks (1 letter per 30-s block) alternated with baseline 
rest blocks (30-s), in which an asterisk (*) was presented for eye fixa- 
tion. Previous studies have demonstrated activation in the anterior 
cingulated and inferior and middle prefrontal cortex in adults (Backes 
et al., 2005; Deblaere et al., 2002; Price, 2010). 
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2.4.2. Motor paradigm 

In the motor task, left handed finger tapping was alternated with 
right handed finger tapping. An arrow was visually presented to indi- 
cate which hand to use. Both tapping conditions were interleaved 
with baseline rest blocks (each block lasts 25 s) and activation maps 
were derived by contrasting the tapping conditions against baseline. 

2.4.3. Functional image processing 

All fMRI time series were processed to generate activation t-maps 
using the Statistical Parametric Mapping (SPM8) (Friston et al, 1994) 
software application (Welcome Department of Cognitive Neurology, 
London, UK), which consisted of several steps. 

First, the images of each dynamic series were realigned to correct for 
head motion using rigid transformations (translation and rotation, 6 de- 
grees of freedom); these motion parameters were compared between 
groups. Next, the Tl -weighted anatomical reference image was segment- 
ed into gray matter, white matter and cerebrospinal fluid components 
and transformed to the Montreal Neurological Institute (MNI) standard- 
ized stereotactic coordinate system. This normalization was subsequently 
applied to the realigned dynamic images. To correct for spatial misregis- 
trations and to strengthen the assumption of normal distribution of 
data, the dynamic images were spatially smoothed using a Gaussian ker- 
nel with a full width at half maximum (FWHM) of 6 mm. 

The task design was convolved with a standard hemodynamic re- 
sponse function to model the fMRI time series. Activation maps were 
statistically evaluated using t-contrasts. 

Second-level analysis was performed to test for differences in over- 
lapping task responses between the two subject groups. 

3. Results 

3.1. Heterogeneity analysis 

The values of the heterogeneity measures are listed in Table 1 for 
both tasks. For the language task, the number of activated regions, frac- 
tal dimension, entropy and homogeneity are slightly higher for patients 
than controls, but the differences were not significant. The isoperimetric 
quotient provided significantly lower values for patients compared to 
controls. Lacunarity showed significantly higher values for patients 
than controls. For the motor task none of the heterogeneity measures 
showed a statistically significant difference between patients and con- 
trols. Fig. 2(a)-(f) shows the heterogeneity measures of patients and 
controls as a function of percentage of most significantly activated 
voxels. All measures showed that the difference between patients and 
controls is conserved and comparable over activation thresholds. 
Fig. 3(a)-(f) shows the heterogeneity measures of patients and controls 
with t-values used as threshold. For increasing a t-threshold less voxels 
are activated. In Figs. 2 and 3 shape and homogeneity values decrease 
with the decreasing numbers of activated voxels, whereas lacunarity, 
entropy and fractals increase. 

The correlations between the different heterogeneity measures are 
listed in Table 2. All measures, other than the number of activation 



regions, show moderate to high correlation (0.7-0.9). Strong negative 
correlation was found between entropy and homogeneity (both based 
on the co-occurrence matrix). 

AUC values for the ROC of each measure are also given in Table 1 . The 
highest AUC was obtained for the lacunarity analysis. The lacunarity 
measure also performed best with a sensitivity of 74% and a specificity 
of 70%. Using the SVM classification, the sensitivity increased to 78% 
and the specificity to 80%. The best two measures were lacunarity 
(F = 0.12) and isoperimetric quotient (F = 0.10). The F-scores of the 
heterogeneity measures are also given in Table 1. The highest 1% 
t-values are concatenated to form a feature vector, which was used for 
the SVM classification. The sensitivity and specificity achieved using 
t-values as features were 45% and 70%, respectively, which further indi- 
cates the importance of using heterogeneity measures. 

3.2. Noise, motion and overlap analysis 

The SNR of the patient group (36.4 ± 0.9) was not significantly 
different from the control group (34.9 ± 1.0). 

Fig. 4(a)-(b) shows the number of clusters and mean cluster size as a 
function of the percentage of the most strongly activated voxels in pa- 
tients and controls. The two figures illustrate that the number of clusters 
is higher in patients than in controls for a range of t-values, whereas 
mean cluster size is lower. To illustrate the distribution of activation 
levels, in Fig. 4(c) the histogram of activation maps for a range of 
t-values is shown. As these distributions do not differ, the total amount 
of activation is comparable between patients and controls. 

Fig. 5(a)-(d) shows the average activation maps of the patient and 
control groups for the language and motor tasks, which were highly 
similar. The group comparison yielded no significant differences. Fig. 6 
shows a group comparison of the Jaccard index between patients and 
controls. The Jaccard index in patients appears to be lower than in con- 
trols indicating that the number of overlapping voxels in patients is 
lower than in controls and the distribution of activated voxels in pa- 
tients is more heterogeneous. 

4. Discussion 

4.1. Current findings 

Different well-known spatial heterogeneity measures were explored 
to capture the complexity, shape and spatial distribution of activated 
brain regions on an individual subject basis. The rationale of the current 
work was to determine spatial heterogeneity measures that can quanti- 
fy the heterogeneity of brain activation. We evaluated the proposed 
measures in patients with epilepsy and language impairment and 
healthy controls in response to a language and motor task. All proposed 
measures of spatial heterogeneity revealed that patients exhibit a more 
heterogeneous activation pattern than controls in response to the lan- 
guage task. Shape- (i.e. isoperimetric quotient) and lacunarity-based 
measures performed better than other methods in distinguishing acti- 
vation patterns of patients from controls. Similar heterogeneity analyses 



Table 1 

Heterogeneity measures for pediatric patients with epilepsy and healthy controls computed with voxels top 5% of t-values, sensitivity and specificity of distinguishing patients and 
controls, the area-under-curve of the receiver-operator characteristic and the F-score. 



Method 




No. of regions 


Fractal dimension 


Isoperimetric quotient 


Entropy 


Homogeneity 


Lacunarity 


SVM 


Language task 


Patients 


44.8 ± 2.9 


2.45 ± 0.001 


0.108 ± 0.0001 


0.280 ± 0.003 


0.973 ± 0.001 


3.22 ± 0.22 






Controls 


38.5 ± 3.6 


2.45 ± 0.001 


0.111 ±0.0001 


0.287 ± 0.002 


0.974 ± 0.001 


2.60 ± 0.13 






p-value 


0.14 


0.79 


0.08 


0.41 


0.50 


0.03 






Sn 


78.2 


86.9 


86.9 


21.7 


34.7 


73.9 


78.3 




Sp 


45.0 


25.0 


45.0 


80.0 


80.0 


70.0 


80.0 




AUC 


0.59 ± 0.01 


0.45 ± 0.01 


0.63 ± 0.02 


0.42 ± 0.01 


0.56 ± 0.01 


0.69 ± 0.02 






F-score 


0.03 


0.01 


0.10 


0.04 


0.01 


0.12 




Motor task 


Patients 


50.8 ± 2.1 


2.90 ± 0.001 


0.0681 ± 0.001 


0.63 ± 0.006 


0.94 ± 0.001 


10.7 ± 0.5 






Controls 


54.5 ± 3.5 


2.90 ± 0.001 


0.0695 ± 0.001 


0.63 ± 0.006 


0.94 ± 0.001 


11.3 ± 0.6 






p-value 


0.42 


0.55 


0.91 


0.60 


0.47 


0.56 
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Fig. 2. Heterogeneity measures in pediatric patients with epilepsy and healthy controls as a function of the percentage of the most strongly activated voxels (3-1 0%). For the fractal di- 
mension (b), lacunarity (f) and entropy (d), the heterogeneity values increase with higher voxel percentages as the complexity, gap size and the distributiveness in the image increase. 
These measures are higher in patients than in controls at all the voxel percentages. The isoperimetric quotient (c) and homogeneity (e) decrease for higher voxel percentages as the ac- 
tivated regions become more irregular and distributed. It can also be observed that over a wide range of thresholds, activated regions in controls always appear closer to the ideal spherical 
shape and more homogenous than in patients. 



were performed for the motor task data. However, for the motor task 
heterogeneity measures did not reveal any differences. This indicates 
that differences in heterogeneity measures are most likely due to the re- 
sponses of the impaired language function. The increased heterogeneity 
was not due to differences in image noise or head motion confounds. 

The Jaccard index in patients was lower than in controls, and the 
spatial distribution of activated voxels in patients is more heteroge- 
neous than in controls. 

The current study provides neuronal correlates, in addition to neuro- 
psychometric measures, for the functional deficits of the epilepsy disorder. 
A goal of such a neuroimaging study is that it adds to the understanding 
how the brain tissue responds differently in children with epilepsy with 
language comorbidity in comparison to healthy controls. As a clinical per- 
spective, this may also serve in future clinical trials as a neuronal biomarker 
to support any treatment and development of the language co-morbidity. 

4.2. Group averaged activation patterns 

The degree of activation map overlap in the patient group was lower 
than in the control group. It was also seen that the number of activated 



regions in patients was higher than in controls, despite the observation 
that standard (overlap based) group comparison revealed no differ- 
ences. These observations suggest that spatial heterogeneity in patients 
is stronger than in controls. However, the quantitative results of the 
overlap and the number of activated regions appeared not significant. 
Moreover, this overlap measure provides a group-based measure, 
which may describe strong variations between subjects of a group but 
not features of spatial heterogeneity in individual subjects. To overcome 
this limitation, we explored measures of spatial heterogeneity that can 
be applied to individual subjects. 

4.3. Nature of spatial heterogeneity in brain activation 

The two most important features in spatial heterogeneity differences 
between individual patients and controls were variations in shape 
(isoperimetric quotient) and sparsity (lacunarity) of activated regions. 
As the isoperimetric quotient was lower in patients, the distribution of 
brain activation regions in patients reflects activated regions with a 
surface area that deviates more strongly from a compact spherical re- 
gion than for healthy controls. As the volume of brain activation was 
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regulated by selecting consistently the highest f-value voxels, the 
shapes of the distributed activated regions are less compact (more 
irregularly shaped and smaller regions) in patients than in controls. 
The observation that the lacunarity measure appears elevated in pa- 
tients demonstrates that activated regions are more separated and 
more irregular than in healthy controls. The number of activated re- 
gions as well as co-occurrence matrix derived measures (i.e. homo- 
geneity and entropy) of spatial heterogeneity were less sensitive to 
differentiate the activation patterns of patients from controls. Co- 
occurrence matrices measure the texture content, in particular the 
repetition of a pattern in an image. The activation patterns in both 
patients and controls are more close to a stochastic texture pattern; 
therefore co-occurrence matrix measures do not provide a strong 
distinction between patients and controls. 

Most remarkable is the observation that activation maps were more 
heterogeneous for the patient group than for the controls for the 
language task. Apparently, the increased heterogeneity of activation 
distribution was related to the impaired language function as for the 
motor task no differences in heterogeneity were found. However, 



more research with more tasks is required to draw final conclusions 
on this matter. 



4.4. Dependence of the different heterogeneity measures 

The combination of measures provides a better discrimination be- 
tween activation patterns of patients and controls than using a single het- 
erogeneity measure. The combination of measures indicates that the 
spatial organization of the activated brain regions is a complex feature 
that cannot be captured completely by one single heterogeneity measure. 
This is also observed in the correlation matrix for all the measures. All 
measures, other than the number of activated regions, were moderately 
correlated because they all capture to some extent the distributiveness 
in the activation map. However, the correlation is not perfect, which 
also points out that one single measure is not sufficient to capture the het- 
erogeneity in patients and controls. This indicates that spatial heterogene- 
ity in fMRI activation maps is due to the variation in all the three 
measures, i.e. shape, gappiness and complexity of the activation regions. 




Fig. 3. Heterogeneity i 
mension (b), lacunar] 
(e) decrease. 



measures in pediatric patients with epilepsy and healthy cons ols with t-values used as threshold on number of voxels activated. Similar to Fig. 2, for the fractal di- 
ity (f) and entropy (d). the heterogeneity value increases with the increase in number of activated voxels and the isoperimetric quotient (c) and homogeneity 
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Table 2 

Table showing the correlation between different measures used in the study for language 
tasks. It can be noted that other than number of activated regions, all the measures show 
moderate to high correlation. 

Fractal Shape Entropy Homogeneity Lacunarity | 1 1 



No. of regions 0.04 -0.15 0.09 -0.05 0.05 



Shape 



Entropy 



Homogeneity 




4.5. Interpretation of increased spatial heterogeneity 

The stronger spatial heterogeneity of activated regions in patients 
compared to controls is observed in shape and spatial distribution of ac- 
tivated regions. This can be explained in the light of brain function. The 
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Fig. 4(a) Number of clusters; (b) mean cluster size for language task as a function of the 
percentage of the most strongly activated voxels in patients and controls and 
(c) histogram of t-values. (a,b) Illustrate that number of clusters are more in patients 
than in controls for a range of t-values whereas mean cluster is less in patients than in con- 
trol. This effect shows that patients' activation pattern is more heterogeneous (distribut- 
ed) than controls, however the effect is not statistically significant using only the 
number of clusters, (c) Shows that statistically the difference in number of activated 
voxels between patients and controls at different intervals of t-values is insignificant. 
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Fig. 5. Activation (overlap) maps of the patient (a,c) and control (b,d) groups for the 
language (a,b) and motor (c,d) tasks. Depicted are the significantly activated voxels in 
color (p < 0.05 ). Similar regions are activated for the language and motor tasks in patients 
and controls. Second level analysis did not provide any statistically significant differences 
between the groups for both tasks. 

brain is organized in functional networks of specialized cortical and 
subcortical regions. One of the interpretations of increased spatial heteroge- 
neity in patients is that auxiliary brain regions may co-activate to compen- 
sate for disturbances in specific functional networks (Dupont et al., 2001 ; 
Eliassen et al, 2008; Vlooswijk et al, 2010). Such a compensation mecha- 
nism creates more irregular and smaller regions in the activation maps 
and also increases spatial dispersion. Another explanation for the increase 
in heterogeneity is that some areas in the impaired brains could activate 
even when that part of the brain is commonly activated for task 
performance (Eliassen et al, 2008). This may affect (suppress) the usual 
part of the brain to activate, and may increase the spatial heterogeneity. 

The selection of a predefined percentage of the strongest activated 
voxels might include more false-positives that are randomly distributed 
throughout the brain. However, analyzing the heterogeneity measures 
by thresholding on activation t-value rather than percentage of activat- 
ed voxels provided similar results. Moreover, the distribution of 
t-values, the image noise and head motion parameters were compara- 
ble for the patient and control groups. 

5. Conclusion 

Novel methods to quantify spatial heterogeneity in brain activation 
maps were explored. These methods focus on different aspects of het- 
erogeneity of the distribution of activated regions, such as shape, 



0.04 



0.03 



■o 0.02 



o.oi 



Patients 

Controls 




2 3 
Threshold on t-values 



Fig. 6. Jaccard index of overlapping voxels in patients and controls as a function of activation 
threshold. 
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complexity, geometric structure and gaps, which can be employed 
at the individual subject level. An increase in spatial heterogeneity 
was observed in pediatric patients with epilepsy relative to healthy 
controls, whereas no significant activation map differences were 
observed when employing conventional (overlap based) group 
comparison. Spatial heterogeneity in patients was best quantified 
in terms of shape and lacunarity measures, describing the distribu- 
tion of the increased surface areas and the irregular gaps of the dis- 
tributed activation pattern, respectively. We propose that spatial 
heterogeneity is a valuable functional (neuroimaging) biomarker 
that could be added to regular brain activation map analysis to char- 
acterize conditions of cerebral disease and should further be ex- 
plored for other neurological disorders as well. Such neuronal 
biomarkers may also help in future trials that aim to diagnose and 
monitor patients that do or not respond to therapeutic assessments 
that aim to improve the language and other cognitive performances. 
We also expect that the applicability of the heterogeneity measures 



can be easily extended to other neuroimaging methods, for instance 
perfusion and diffusion imaging. 
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Appendix A 

Conceptual set of images (dimension, 288 x 288) to illustrate vari- 
ous aspects of spatial heterogeneity in 2D are shown in Fig. A.l. From 
image (a) to (f) the spatial heterogeneity increases in terms of number 
of regions, organization of regions and variation in region size and 




Image 


(a) 


(b) 


(c) 


(d) 


(e) 


(f) 


No. of regions 


1 


9 


9 


9 


9 


5582 


Fractal 
dimension 


1.98 


2.59 


2.65 


2.41 


2.46 


2.45 


Shape 


0.70 


0.08 


0.08 


0.10 


0.09 


0.00094 


Entropy 


1.31 


1.70 


1.75 


1.60 


1.69 


4.01 


Homogeneity 


3094 


2899 


2829 


2937 


2906 


1661 


Lacunarity 


2.11 


9.94 


11.40 


11.96 


18.96 


33 



Fig. A.l. Conceptual set of images (dimension, 288 x 288) to illustrate various aspects of spatial heterogeneity in 2D. From image (a) to (f) the spatial heterogeneity increases in terms of 
number of regions, organization of regions and variation in region size and shape. The set of model images in combination with the table demonstrates the change in heterogeneity mea- 
sures with the various types of region distributions. 
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shape. Image (f) shows one circular activated region and number of small 
randomly activated regions scattered around the complete image, which 
represent noise. The set of model images in combination with the table 
demonstrates the change in heterogeneity measures with the various 
types of region distributions. Each image contains the same number of 
white (activated) pixels, thus the percentage of activated voxels is con- 
stant over images. 

The lacunarity (i.e. gappiness) increases from image (a) to (f) as the 
gaps between the regions increase (in size and variation) from (a) to (f). 
The grid size used to compute lacunarity is 1 6 x 1 6 pixels. 

Fractal dimension, which is a measure of complexity, increases from 
(a) to (c) because the complexity of the image increases due to modifi- 
cations in shape and position of the regions. From image (c) to (d), a few 
regions have increased in size, which has reduced the complexity in the 
image thus reducing the fractal dimension compared to image (c). Frac- 
tal dimension in image (e) is higher than in (d), because the complexity 
of the objects has increased, however it is less than in image (c) because 
of the few larger regions in image (e) compared to (c). 

Shape (in terms of the isoperimetric quotient) of the regions in image 
(a) is closest to a perfect sphere, and therefore the shape measure is rela- 
tively high (towards the value of one). By dividing the large circular re- 
gion of image (a) into smaller regions, the perimeter of the regions 
increases and the shape measure decreases. The shape of the regions re- 
mains constant from image (b) to (c). It further increases in image 

(d) as a few regions have become larger than in image (c), which reduces 
the perimeter. It is smaller in image (e) compared to (d) because the re- 
gions are now elliptical in shape. It is the smallest in (f) because of 
many small regions in the image. 

Entropy is the measure of randomness, thus random distribution of 
activated voxels. It is the lowest in image (a) because of the homogenous 
surface. It shows the higher value in images (c) and (d) because all the 
(activated) regions are randomly distributed thus creating randomness 
in the image. Homogeneity is opposite to entropy. It shows the higher 
value in image (a) compared to images in (c) and (d). Entropy is maximal 
and homogeneity is minimal in (f) because of the extensive noise. 

The number of (activated) regions is the number of disconnected 
components or regions in an image. They are two in (a), nine in (b) to 

(e) and 5582 in (f). 

Appendix B 

6. Motion analysis: Displacement, smoothness and entropy 

6.1. Displacement 

The average Euclidean (i.e. root-mean-squared) distance for the two 
groups is shown in Table B.l. The difference between two groups is not 
significant. 

6.2. Smoothness 

The smoothness of the displacement parameter time-series was 
expressed in terms of its bandwidth. Smoothness results are provided 
in Table B.l for the two groups. Statistical analysis revealed no differ- 
ences between patients and controls. 

6.3. Entropy 

Entropy is a measure of randomness in a signal. Entropy of the dis- 
placement parameter of the two groups is listed in Table B.l. Statistically 
there is no significant difference between the two groups. 

6.4. Correlation analysis 

Displacement, smoothness and entropy poorly correlate (< 0.3 ) with 
all the heterogeneity parameters proposed, see Table B.2. 

These analyses show that head motion was not different between 
patients and controls and does not relate to the heterogeneity of the ac- 
tivation maps. 
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Table B.l 

Displacement, smoothness and entropy parameters of the patient and control groups for 
the language task. None of the comparisons appeared statistically significant. 





Displacement 


Smoothness 


Entropy 






Absolute Relative 
(mm) (mm) 


Absolute Relative 

(s- 1 ) (s->) 


Absolute 


Relative 


Patient 
Control 


0.57 ± 0.34 0.13 ±0.16 
0.65 ±0.38 0.13 ±0.17 


3.2 ±2.2 6.4 ±3.6 
2.9 ± 1.3 6.9 ±3.8 


7.32 ± 0.12 
7.42 ± 0.10 


6.99 ± 0.45 
7.12 ± 0.27 



Table B.2 

Correlation between heterogeneity measures and displacement, smoothness and entropy. 
All heterogeneity measures poorly correlate with displacement, smoothness and entropy. 
These measures are computed on the language task. 



Displacement Smoothness Entropy 





Absolute 


Relative 


Absolute 


Relative 


Absolute 


Relative 


Number of regions 


0.27 


0.31 


0.32 


0.21 


-0.20 


-0.31 


Fractal dimension 


0.12 


0.16 


0.1 1 


-0.09 


-0.17 


-0.18 


Shape 


-0.04 


-0.18 


-0.02 


0.19 


0.19 


0.14 


Entropy 


0.10 


0.15 


0.09 


-0.02 


-0.24 


-0.23 


Homogeneity 


-0.08 


-0.10 


-0.07 


0.01 


0.22 


0.22 


Lacunarity 


0.01 


0.3 


0.02 


-0.26 


-0.17 


-0.04 
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